Motivation

Climate issues are pervasive but typically disproportionately affect low income communities and developing countries. Our group wanted to explore how air pollution has changed over time and affect countries differently. Specifically, we wanted to analyze how a country’s economic and social position can either increase, decrease, or not have observable impact on the affects of air pollution. In laymen terms, does air pollution affect underdeveloped countries disproportionately?

Set Up

Before we start, we need to ensure that we have all the relevant libraries installed and imported.

Run these in the console, or only the ones that your system does not have, to install packages in addition to the ezids package.

install.packages("tidyverse")
install.packages("rworldmap")
install.packages("tmap")
install.packages("spData")
install.packages("sf")
install.packages("ggpubr")
install.packages("dplyr")
install.packages("knitr")
install.packages("magrittr")

Data Sources

For our analysis, we will be working with 5 main data sources shown in the table below:

The main variables in our datasets will include:

Data Wrangling

While data from Kaggle are already in a format to be cleaned, downloaded data from United Nations required a little data wrangling. Mainly, we needed to extract just countries’ data from the Excel workbooks and into their own contained csv files. Since we only need to do this once and programming it would take significant time to choose the specific cells that we need, we opted to perform this step outside of R and in Excel. Note that if this were a part of a real production data pipeline, we would take the time to program the data extraction but would likely choose a different programming language such as Python that is a bit more robust in these types of tasks like web scraping and data transformations in Pandas.

Figure 3: Sample screenshot of data downloaded from UN including unnecessary elements like banners and other regional data.

Figure 4: Sample screenshot of transformed UN dataset.

Load Data

Clean Data

First thing that we need to drop unnecessary columns and set datatypes (factor, num, etc.).

Do lower GDP countries have more deaths per 100k due to air pollution?

## 
## Call:
## lm(formula = log(Deaths.Air.Pollution.per.100k) ~ log(gdp.per.capita), 
##     data = final_df_sf)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.099 -0.235  0.000  0.206  1.431 
## 
## Coefficients:
##                     Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)          7.38778    0.02663     277 <0.0000000000000002 ***
## log(gdp.per.capita) -0.38952    0.00323    -121 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.369 on 5195 degrees of freedom
## Multiple R-squared:  0.737,  Adjusted R-squared:  0.737 
## F-statistic: 1.45e+04 on 1 and 5195 DF,  p-value: <0.0000000000000002

## $`1`
## [1]  95.2 920.3
## 
## $`2`
## [1]  921 3100
## 
## $`3`
## [1]  3104 11493
## 
## $`4`
## [1]  11497 119106

##               Df   Sum Sq Mean Sq F value              Pr(>F)    
## qnt            3 10735714 3578571    3133 <0.0000000000000002 ***
## Residuals   5193  5932162    1142                                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1. What are the impacts of population size from GDP and Deaths due to Air Pollution globally?

1. Categorizing GDP into low (0), medium(1) and high(2)

SQ6 <- final_df
library(dplyr)
groupings_pop <- SQ6 %>% mutate(population_grouping = case_when(Population.thousands >=  100001 & Population.thousands <= 1000000000 ~ 2,
                                             Population.thousands >= 50001  & Population.thousands <= 100001 ~ 1,
                                             Population.thousands >= 0 & Population.thousands <= 50000 ~ 0)) # end function

2. Building a logit model GDP as predictor on Population

pop_gdpLogit <- glm(population_grouping ~ groupings_pop$gdp.per.capita, data = groupings_pop)

summary(pop_gdpLogit)

## 
## Call:
## glm(formula = population_grouping ~ groupings_pop$gdp.per.capita, 
##     data = groupings_pop)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -0.246  -0.185  -0.181  -0.180   1.820  
## 
## Coefficients:
##                                 Estimate  Std. Error t value
## (Intercept)                  0.180177461 0.008492388   21.22
## groupings_pop$gdp.per.capita 0.000000553 0.000000456    1.21
##                                         Pr(>|t|)    
## (Intercept)                  <0.0000000000000002 ***
## groupings_pop$gdp.per.capita                0.23    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.267)
## 
##     Null deviance: 1389.8  on 5196  degrees of freedom
## Residual deviance: 1389.4  on 5195  degrees of freedom
## AIC: 7899
## 
## Number of Fisher Scoring iterations: 2

3. Bulding a logit model GDP and Deaths due to Air Population as predictor on Population

pop_gdp_deathsLogit <- glm(population_grouping ~ groupings_pop$gdp.per.capita + groupings_pop$Deaths.Air.Pollution.per.100k, data = groupings_pop)

summary(pop_gdp_deathsLogit)

## 
## Call:
## glm(formula = population_grouping ~ groupings_pop$gdp.per.capita + 
##     groupings_pop$Deaths.Air.Pollution.per.100k, data = groupings_pop)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -0.215  -0.194  -0.186  -0.169   1.842  
## 
## Coefficients:
##                                                 Estimate   Std. Error t value
## (Intercept)                                  0.202910224  0.018188095   11.16
## groupings_pop$gdp.per.capita                 0.000000136  0.000000543    0.25
## groupings_pop$Deaths.Air.Pollution.per.100k -0.000213150  0.000150811   -1.41
##                                                        Pr(>|t|)    
## (Intercept)                                 <0.0000000000000002 ***
## groupings_pop$gdp.per.capita                               0.80    
## groupings_pop$Deaths.Air.Pollution.per.100k                0.16    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.267)
## 
##     Null deviance: 1389.8  on 5196  degrees of freedom
## Residual deviance: 1388.9  on 5194  degrees of freedom
## AIC: 7899
## 
## Number of Fisher Scoring iterations: 2

2. What are the effects of (low or high) GDP and population size on deaths due to air pollution in Sub-Saharan Africa?

3. Can we let the data tell us what type of groupings exist in our dataset? How consistent are they to our preconceived groupings such as region or developed vs developing countries?

We wanted to better understand what types of groupings exist within our dataset. More concretely, we wanted to put aside our preconceived presumptions about our dataset and create a K-Means clustering model and allow the the data to guide us to the possible clusters that may exist within our data.

K-Means clustering is an unsupervised machine learning algorithm that is very useful to parse the dataset and identify K number of groups where observations within each group have high similarity.

First, let’s create an index label.

Next, let’s select only the numerical values from our dataset.

4. Can we make a prediction of the future GDP by considering the population size, location, and the reation of population and deaths due to air pollution?

For our major model we make a nice prediction for GDP through total deaths caused by air pollution, so as we can see that through the distribution of the data, there are strong relationships between total deaths caused by air pollution, deaths caused by indoor air pollution, and deaths caused by outdoor air pollution.

# select what we need
# total
#ggplot(total3, aes(x=Year, y=air_pollution_total)) +
#    geom_bar(stat='identity', position='dodge')

# indoor
#ggplot(total3, aes(x=Year, y=air_pollution_indoor)) +
#    geom_bar(stat='identity', position='dodge')
# outdoor particulate
#ggplot(total3, aes(x=Year, y=air_pollution_outdoorP)) +
#    geom_bar(stat='identity', position='dodge')
# outdoor ozone 
#ggplot(total3, aes(x=Year, y=air_pollution_outdoorO)) +
#    geom_bar(stat='identity', position='dodge')


air_pollution_df_cleaned1 <- air_pollution_df

So it seems like we can also get a nice GDP prediction by using deaths caused by indoor air pollution and deaths caused by outdoor air pollution as the basis of the prediction model.

# make Entity and year factor
air_pollution_df_cleaned1$Entity = factor(air_pollution_df_cleaned1$Entity)
air_pollution_df_cleaned1$Code = factor(air_pollution_df_cleaned1$Code)
air_pollution_df_cleaned1$Year = factor(air_pollution_df_cleaned1$Year)

# polution
air_pollution_df_cleaned1 <- rename(air_pollution_df_cleaned1, Country = Entity, ISO.alpha3.code = Code, Deaths.Air.Pollution.per.100k = Air.pollution..total...deaths.per.100.000. , Deaths.Air.Pollution.Indoor.per.100k = Indoor.air.pollution..deaths.per.100.000., Deaths.Air.Pollution.OutdoorP.per.100k = Outdoor.particulate.matter..deaths.per.100.000., Deaths.Air.Pollution.OutdoorO.per.100k = Outdoor.ozone.pollution..deaths.per.100.000.)


# population
#population_region_df_cleaned <- subset(population_region_df, select = -c(Notes, Type, Parent.code))

#population_region_df_cleaned <- population_region_df_cleaned %>%
#  pivot_longer(
#    cols = starts_with("X"), 
#    names_to = "Year", 
 #   values_to = "Population.thousands", 
  #  values_drop_na = TRUE
#  )

# population_region_df_cleaned$Year<-gsub("X","",as.character(population_region_df_cleaned$Year))
# population_region_df_cleaned <- as.data.frame(apply(population_region_df_cleaned, 2, function(x) gsub("\\s+", "", x))) 

# population_region_df_cleaned$Country = factor(population_region_df_cleaned$Country)
# population_region_df_cleaned$SDGRegion = factor(population_region_df_cleaned$嚜燙DGRegion)
# population_region_df_cleaned$Country.code = factor(population_region_df_cleaned$Country.code)
# population_region_df_cleaned$SubRegion = factor(population_region_df_cleaned$SubRegion)
# population_region_df_cleaned$Year = factor(population_region_df_cleaned$Year)
# population_region_df_cleaned$Population.thousands = as.numeric(population_region_df_cleaned$Population.thousands)

# population_region_df_cleaned <- rename(population_region_df_cleaned, M49.code = Country.code)

# country_codes_df$M49.code = factor(country_codes_df$M49.code)
# country_codes_df$Country.or.Area = factor(country_codes_df$嚜澧ountry.or.Area)
# country_codes_df$ISO.alpha3.code = factor(country_codes_df$ISO.alpha3.code)
# country_codes_df$ISO.alpha2.code = factor(country_codes_df$ISO.alpha2.code)

# world$iso_a2 = factor(world$iso_a2)
# 
# world_df_cleaned <- subset(world, select = c(iso_a2, geom))

# merge
final_dfc <- merge(x = air_pollution_df_cleaned1, y = country_codes_df, by = 'ISO.alpha3.code')

final_dfc <- merge(x = final_dfc, y = gdp_df_cleaned, by = c('ISO.alpha3.code', 'Year'))
final_dfc <- merge(x = final_dfc, y = population_region_df_cleaned, by = c('M49.code', 'Year'))
final_dfc <- merge(x = final_dfc, y = world_df_cleaned, by.x = 'ISO.alpha2.code', by.y = 'iso_a2', all.x = TRUE) #left join

# Remove columns
final_dfc <- subset(final_dfc, select = -c(Country.or.Area, Country.y, Country))

# Calculate GDP per capita 
final_dfc$gdp.per.capita <- final_dfc$GDP.USD / final_dfc$Population.thousands


final_dfc$Deaths.Air.Pollution.Outdoor.Total.per.100k <- final_dfc$Deaths.Air.Pollution.OutdoorO.per.100k + final_dfc$Deaths.Air.Pollution.OutdoorP.per.100k

For the first model, we build a linear model using ‘gdp.per.capita’ and ‘Deaths.Air.Pollution.Indoor.per.100k’.

final_df_sfc = st_as_sf(final_dfc)
fit2c <- lm(gdp.per.capita ~Deaths.Air.Pollution.Indoor.per.100k, data=final_df_sfc)

summary(fit2c)

## 
## Call:
## lm(formula = gdp.per.capita ~ Deaths.Air.Pollution.Indoor.per.100k, 
##     data = final_df_sfc)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -15508268  -8789754  -3652222   3500260 102785447 
## 
## Coefficients:
##                                      Estimate Std. Error t value
## (Intercept)                          16339920     254681    64.2
## Deaths.Air.Pollution.Indoor.per.100k  -126647       3308   -38.3
##                                                 Pr(>|t|)    
## (Intercept)                          <0.0000000000000002 ***
## Deaths.Air.Pollution.Indoor.per.100k <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13900000 on 5195 degrees of freedom
## Multiple R-squared:  0.22,   Adjusted R-squared:  0.22 
## F-statistic: 1.47e+03 on 1 and 5195 DF,  p-value: <0.0000000000000002

ggplot(final_df_sfc, aes(gdp.per.capita, Deaths.Air.Pollution.Indoor.per.100k)) +
  geom_point() +
  geom_smooth(method='lm', se=FALSE) +
  stat_regline_equation(label.y = 350, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 400, aes(label = ..rr.label..))

Figure 6.4-1: Deaths due to Indoor Air Pollution vs GDP per capita.

The R-squared value of this model is 0.22, which is a really bad result. So follow the steps in the major model building, now we put these two features into a log() and it should performed better.

## 
## Call:
## lm(formula = log(gdp.per.capita) ~ log(Deaths.Air.Pollution.Indoor.per.100k), 
##     data = final_df_sfc)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.354 -0.599  0.098  0.647  2.892 
## 
## Coefficients:
##                                           Estimate Std. Error t value
## (Intercept)                               16.33257    0.01763     927
## log(Deaths.Air.Pollution.Indoor.per.100k) -0.54761    0.00517    -106
##                                                      Pr(>|t|)    
## (Intercept)                               <0.0000000000000002 ***
## log(Deaths.Air.Pollution.Indoor.per.100k) <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.891 on 5195 degrees of freedom
## Multiple R-squared:  0.684,  Adjusted R-squared:  0.684 
## F-statistic: 1.12e+04 on 1 and 5195 DF,  p-value: <0.0000000000000002

Figure 6.4-2: Log Deaths due to Indoor Air Pollution vs Log GDP per capita.

As we can see, although it performed much better with a R-squared value 0.68 but this is still not good as the major model.

Next we try the deaths caused by outdoor air pollution ‘Deaths.Air.Pollution.Outdoor.Total.per.100k’.

## 
## Call:
## lm(formula = gdp.per.capita ~ Deaths.Air.Pollution.Outdoor.Total.per.100k, 
##     data = final_df_sfc)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -14051236  -9537840  -5335898   2518144 106063637 
## 
## Coefficients:
##                                             Estimate Std. Error t value
## (Intercept)                                 15643215     461880    33.9
## Deaths.Air.Pollution.Outdoor.Total.per.100k  -150382      10830   -13.9
##                                                        Pr(>|t|)    
## (Intercept)                                 <0.0000000000000002 ***
## Deaths.Air.Pollution.Outdoor.Total.per.100k <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15400000 on 5195 degrees of freedom
## Multiple R-squared:  0.0358, Adjusted R-squared:  0.0356 
## F-statistic:  193 on 1 and 5195 DF,  p-value: <0.0000000000000002

Figure 6.4-3: Deaths due to Outdoor Air Pollution vs GDP per capita.

The performence of this model is also bad, with a R-squared value 0.036, so we put the features in the log() too.

## 
## Call:
## lm(formula = log(gdp.per.capita) ~ log(Deaths.Air.Pollution.Outdoor.Total.per.100k), 
##     data = final_df_sfc)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.591 -1.272 -0.027  1.297  3.471 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                       15.6988     0.1565  100.28
## log(Deaths.Air.Pollution.Outdoor.Total.per.100k)  -0.1991     0.0442   -4.51
##                                                              Pr(>|t|)    
## (Intercept)                                      < 0.0000000000000002 ***
## log(Deaths.Air.Pollution.Outdoor.Total.per.100k)            0.0000068 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.58 on 5195 degrees of freedom
## Multiple R-squared:  0.00389,    Adjusted R-squared:  0.0037 
## F-statistic: 20.3 on 1 and 5195 DF,  p-value: 0.00000675

Figure 6.4-4: Log Deaths due to Outdoor Air Pollution vs Log GDP per capita.

A surprised point is that although we put the variables into log(), the performance of model is still bad and even worse. The R-squared value of this model is 0.0039.

Finally, we are going to test build a model base on both deaths caused by outdoor and indoor air pollution. To check if this model also works bad.

## 
## Call:
## lm(formula = gdp.per.capita ~ Deaths.Air.Pollution.Indoor.per.100k + 
##     Deaths.Air.Pollution.Outdoor.Total.per.100k, data = final_df_sfc)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -18198746  -7911990  -2886359   3411956  96919910 
## 
## Coefficients:
##                                             Estimate Std. Error t value
## (Intercept)                                 26403084     457576    57.7
## Deaths.Air.Pollution.Indoor.per.100k         -144489       3190   -45.3
## Deaths.Air.Pollution.Outdoor.Total.per.100k  -242555       9394   -25.8
##                                                        Pr(>|t|)    
## (Intercept)                                 <0.0000000000000002 ***
## Deaths.Air.Pollution.Indoor.per.100k        <0.0000000000000002 ***
## Deaths.Air.Pollution.Outdoor.Total.per.100k <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13100000 on 5194 degrees of freedom
## Multiple R-squared:  0.309,  Adjusted R-squared:  0.309 
## F-statistic: 1.16e+03 on 2 and 5194 DF,  p-value: <0.0000000000000002

Figure 6.4-5: Deaths due to Indoor + Outdoor Air Pollution vs GDP per capita.

The R-squared value of this model is 0.29. Next we put these variables into log().

## 
## Call:
## lm(formula = log(gdp.per.capita) ~ log(Deaths.Air.Pollution.Outdoor.Total.per.100k) + 
##     log(Deaths.Air.Pollution.Indoor.per.100k), data = final_df_sfc)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.200 -0.594  0.055  0.621  2.901 
## 
## Coefficients:
##                                                  Estimate Std. Error t value
## (Intercept)                                      17.60605    0.08820   199.6
## log(Deaths.Air.Pollution.Outdoor.Total.per.100k) -0.35985    0.02444   -14.7
## log(Deaths.Air.Pollution.Indoor.per.100k)        -0.55212    0.00507  -108.8
##                                                             Pr(>|t|)    
## (Intercept)                                      <0.0000000000000002 ***
## log(Deaths.Air.Pollution.Outdoor.Total.per.100k) <0.0000000000000002 ***
## log(Deaths.Air.Pollution.Indoor.per.100k)        <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.873 on 5194 degrees of freedom
## Multiple R-squared:  0.696,  Adjusted R-squared:  0.696 
## F-statistic: 5.96e+03 on 2 and 5194 DF,  p-value: <0.0000000000000002

Figure 6.4-6: Log Deaths due to Indoor + Outdoor Air Pollution vs Log GDP per capita.

For this model the R-squared value is 0.69. This is actually not a bad result but when comparing to the major model using total deaths caused by air pollution, 0.69 is not good enough for predicting GDP.

We can get the conclusion that we can not get a nice GDP prediction by using deaths caused by indoor air pollution and deaths caused by outdoor air pollution as the basis of the prediction model.